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ABSTRACT 

We describe the first results from a six-month long reverberation-mapping experiment in the ultravi¬ 
olet based on 171 observations of the Seyfert 1 galaxy NGC 5548 with the Cosmic Origins Spectrograph 
on the Hubble Space Telescope. Significant correlated variability is found in the continuum and broad 
emission lines, with amplitudes ranging from ^ 30% to a factor of two in the emission lines and a 
factor of three in the continuum. The variations of all the strong emission lines lag behind those of the 
continuum, with HeiiA1640 lagging behind the continuum by ^ 2.5 days and LyaA1215, CivA1550, 
and SiivA1400 lagging by ^ 5-6 days. The relationship between the continuum and emission lines 
is complex. In particular, during the second half of the campaign, all emission-line lags increased by 
a factor of 1.3-2 and differences appear in the detailed structure of the continuum and emission-line 
light curves. Velocity-resolved cross-correlation analysis shows coherent structure in lag versus line-of- 
sight velocity for the emission lines; the high-velocity wings of Civ respond to continuum variations 
more rapidly than the line core, probably indicating higher velocity BLR clouds at smaller distances 
from the central engine. The velocity-dependent response of Lya, however, is more complex and will 
require further analysis. 

Subject headings: galaxies: active — galaxies: individual (NGC 5548) — galaxies: nuclei — galaxies: 
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1. INTRODUCTION 


I.l. The Broad-Line Region 


One of the most prominent characteristics of the ul¬ 
traviolet (UV), optical, and near-infrared (NIR) spec¬ 
tra of active galactic nuclei (AGNs) is the presence of 
broad emission lines. While we know that these fea¬ 
tures arise on scales not much larger than the accretion 
disk, their physical nature remains one of the major un¬ 
solved mysteries in AGN astrophysics. A particularly 
important feature of the broad emission lines is that they 
are, by definition, resolved in line-of-sight (LOS) veloc¬ 
ity, and their large widths leave little doubt that the pri¬ 
mary broadening mechanism is differential Doppler shifts 
due to the motion of individual gas clouds, filaments, or 
more-or-less continuous flows around the central black 
hole. 

However, it is not possible to establish the broad¬ 
line region (BLR) kinematics simply by inverting the 
line profiles because this inverse problem is degener¬ 
ate, with a wide variety of si mple velocity models pro¬ 
vidin g satisfactory fits (e.g., ICaoriotti. Foltz. &: BvardI 
Il980f) . The existing evidence on the BLR kinemat¬ 
ics is ambiguous: some of this gas may flow in¬ 
ward, helping to feed the central black hole. Ex¬ 
tended, flattened, rotating disk-like structures seem to 
be important in at least some BLRs, as shown sta¬ 
tistically for radio-loud AGNs (I Wills fc Brownd 119861: 
Vestergaard. Wilkes, fc Barthell 120001: IJ arvis fc McLurel 
200011 . by the pronounced doubl e-peaked profiles ob¬ 
serve d in s ome sources /e.g.. lEracleou s fc Hal oernl 
1994 I 2 OO 3 I: iStrateva et al.l 120031 iGezari et al.1 l2007t 
Lewis. Eracleous. fc Storc hi-B ergmann 1201011 , ^d from 
spect ropolarimetry (|Smith et al.l 12004 IVoung et al.l 
l2007fl . There is evidence of the importance of the black 
hole gravity in domin ating the motion of the BLR gas 
(|Peterson et al.l 120041. although rad i ation pressure may 
also p lay a role ((Marconi et al.l[2008t iNetzer fc MarziimJ 

Imof) . 

On the other hand, developments over the last 
two decades re-open the interesting possibility that 
much of the e mitting BLR gas i s due to outflow¬ 
ing w i nds fe.g.. iBottorff et al.lll997t iMnrrav fc Ghiangl 
IProga etaLI 120001: lEverettI 120031 IEIvIfJ 12004 
Young et al.l l200’?fl . perhaps connecte d to the outflows 
detec t ed in absorption featu re s (e.g .. iHamann fc Sabrd 
2 OOI ll^gold et all [2?M I2(M iKriss et al.l l2?)Tir 
Kaastra et al.ll20l4lScott et al.ll2014 , whose kinematics 
and energetics are also poorly understood. The unknown 
dynamics of the BLR gas represents a serious gap in our 
understanding of AGNs and in the calibrations needed 
for the study of black-hole/host-galaxy co-evolution up 
to very high redshifts. 

There have been many attempts to model the physics 
of the BLR. In general, photoionization equilibrium mod¬ 
els can reproduce the line intensities, but self-consistent 
models that provide simultaneous solutions to the line in¬ 
tensities, profiles, and variability ar e lacking. The l ocally 
optimally emitt i ng cl oud model ([Baldwin et al.l 119951 : 
iKorista fc GoadI 2000|fl and the stratified cloud model 
((Kaspi fc Netzeij 1999( 1 explain most observed line inten¬ 
sities and some of the observed time lags between the 
continuum and emission lines. However, they lack the 
important kinematic ingredients required to explain the 


observed line profiles. 


1.2. Reverberation Mapping 

In order to understand the structure and kinematics 
of the BLR, we must break the degeneracy that comes 
from the study of the line profiles alone. We can do 
this by using reverberation mapping (RM) to determine 
how gas at various LOS velocities responds to contin- 
uum variations as a funct i on of light travel-tim e delay 
((Blandford fc McKeelll98d: lPeterso^ll993l [20141 . 

Over the last quarter century, the RM technique has 
become a standard tool for investigating the BLR. In 
its simplest form, RM is used to determine the mean 
time delay between continuum and emission-line varia¬ 
tions, typically by cross-correlation of the respective light 
curves. It is assumed that this represents the mean light- 
travel time across the BLR. By combining this with the 
emission-line width, which is assumed to reflect the ve¬ 
locity dispersion of gas whose motions are dominated by 
the mass of the central black hole, the black hole mass 
can be estimated. RM in this form has been used to mea¬ 
sure the black h ole masses in over 50 AGNs (for a recent 
compilation, see iBentz fc Katd 12015( 1 to a typical accu¬ 
racy of 0.3 dex. Important findings that have arisen 
from these RM studies include the following: 

1. In a given AGN, emission lines that are character¬ 
istic of higher-ionization gas respond more rapidly 
to continuum flux variations than those character¬ 
istic of lower-ionization gas, in dicating ionization 
stratification within the BLR ((Glavel et al.l 119911 : 
l^ichert et ^11994 1. 


2. There is an inverse correlation between the time 
delay, or lag r, for a particular emission line and 
the Doppler width AU of that emission line. The 
relationship for a given AGN is consistent with the 
virial prediction AU oc jlPeterson_&_Wand^ 
I 999 LI 2 OOOI: ^llatschnvll2003l: IPeterson et al.ll2004 


Bentz et al. I2010all . Without this relationship, RM 

masses would be highly dubious. 


3. There is an empirical relationship between the 
AGN luminosity L and the radius of the BLR 
R (hereafter the R-L relationship) that is 
well-established only fo r the H/3 emission line 
(iKaspi et al.l [200(1 [20M IBentz et al.l [20M [^(M 
I 2 OI 3 II . Limited data on GivA1549 indicates 
a si milar relationship applies to that line as 
well (iPeters on et al.l l2005l: (Ves tergaar d fc PetersonI 

[20?)5lKasDi et al.ll2nn7l : iPark et al.ll2ni3ll . The ex- 

istence of R-L relationships for both low-ionization 
and high-ionization lines has been independently 
confir med by gravitationa l microlensing observa¬ 
tions (iGuerras et al.l[20T^ . 

The R-L relationship is of particular interest as it allows 
estimation of the central black-hole mass based on a 
single spectrum from which the line width is measured 
and the BLR radius is inferred from the AGN luminosity. 
This neatly bypasses the need for a direct RM measure¬ 
ment of the emission-line time lag. RM is necessarily 
resource intensive: even to determine the mean time 
delay for an emission line typically requires some 30-50 
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well-spaced high-quality spectrophotometric observa¬ 
tions or a good measure of luck for fewer observations. 
The R-L relationship is very important as the RM-based 
mass deter minations anchor empiri c al scaling re l ation- 


ships fe.ff.. IMcLiire JarvisI 120021 IVesterEraardI 120021: 

Shields et al.l 20031: IGriioe fc MathiiH 20041: Vestersaard 

2004 Greene & Hoi 120051: IMathur & Gruod 

2005t 

Kollmeier et al. 120061: Vestergaard fc Petersonl 

2 ^ 

Salviander et al.l 120071: iTreu et al. 2007: IMcGill et al. 

2008 

: IPark et al.l 120131 120151: INetzer fc Trakhtenbrot 

2014 

) that are used to estimate the masses of quasar 


black holes in large numbers (e.g ., iVestergaard et all 
2008t iVestergaard fc Osmed 120091 : iShen et al.l 1201 iF 
De Rosa et al. I2014D. Virtually all quasar mass esti¬ 


mates and their astrophysical uses are tied to RM. 

Measurement of the mean lag and line width for a given 
emission line provides important, though limited, infor¬ 
mation about the BLR and the central mass of the AGN. 
We are only now beginning to realize the full power of 
RM through velocity-resolved investigations of the BLR 
response. The first generation of successful RM programs 
provided sufficient understanding of AGN variability and 
BLR response times to design programs that could effec¬ 
tively extract velocity-dependent information that would 
lead to an understanding of the structure and kinematics 
of the BLR thro ugh recovery of “v elocity-delay” maps 
from RM data (iHorne et al.l l2004f) . The relationship 
between the continuum variations AC{t) and velocity- 
resolved emission-line variations AL{V,t) is usually de¬ 
scribed as 

nOO 

AL{V,t)= 4'(V, r)AC'(t - r)dr, (1) 
Jo 

where Tfld r') is the “respons e function,” or velocity- 
delay map (jHorne et al.ll200'i l. As can be seen by in¬ 
spection, 4'(y,r) is simply the observed emission-line 
response to a delta-function continuum outburst. The 
velocity-delay map is simply the BLR geometry and 
kinematics projected into the two observable quantities 
of LOS velocity and time delay relative to the contin¬ 
uum. This linearized echo model is justified by the fact 
that the continuum and emission-line variations are gen¬ 
erally q uite small (10-20%) on r everberation time scales 
(see also lGackett fc Hornell2006[) . The technical goal of a 
reverberation program such as the one described here is 
to recover the velocity-delay map 4/(1/, r) from the data 
and thus infer the geometry and kinematics of the BLR. 

Time-resolved velocity-delay maps have now been ob- 
tained for a handful of AGNs (e.g., [Bentz et al.ll2010bl: 
Brewer_et_al. 120111: JPancoast et alJ l2012t iGrier et alJ 
201 ,‘lHPancoast et al.ll2014li . but only for optical lines (the 
Balmer lines, HeiA5876, and HeiiA4686). In general, 
these suggest flattened geometries at small to modest in¬ 
clinations and some combination of virialized motion and 
infall. An outflow si gnature has been o bserved in only 
one case, NGG3227 (iDennev et al.ll200^ . 

The lack of velocity-delay maps for UV lines, on the 
other hand, leaves us with a very incomplete understand¬ 
ing of the BLR. It is, in fact, the high-ionization level UV 
resonance lines (e.g., GivA1549, SiivA1400, LyaA1215) 
that might be expected to dominate any outflowing com¬ 
ponent of the BLR. The optical lines, in contrast, gener¬ 
ally seem to arise in disk-like structures with infall com¬ 


ponents le.g.. iPancoast et 311120141) . 

RM studies in the UV have been limited. Several 
observing campaigns were undertaken with the Inter¬ 
national Ultraviolet Explorer (lUE) or Hubble Space 
Telescope (HST) or both on (i) NGG554 8 (iGlavel et al ' 
1991 : lKorista et aT1ll995h . ( ii) NGG3783 (iReichert et~ 
1^94), (iii) Fairall9 (ICl avel . Wamsteker. fc Glas; 


1989: iRodrfauez-Pascual et al.ni997D . (iv) 3C 390.3 


(lO’Brien et al.l I1998D . (v) NGC 7469 (iWanders et al. 


1997) . ('vi'l NGG4151 (Iciavel et al.lll99fll : Ullrich fc Home 


1^96 : icVenshaw et al.lll996ll . ^ (vij) Akn564 jcollier et~. 

2001 ). and (viii) NGG4395 (jPeterson et al.ll2005D . With 
the exception of Akn564, which showed essentially no 
emission-line variability over a comparatively short 
campaign, all of these programs yielded emission¬ 
line lags, but only limited information about the 
detailed response of th e UV e mission lines (e.g. , 
Horne. Welsh, fc Petersonl 119911: Kroli k et al.l 119911 : 
Wanders et al. l ll995l:lDone fc Kroliklll996D . The existing 
velocity-delay map for NGG4151 shows some incipient 
structure in CivA1549 and HeiiA1640 and a general 
shape that seems to b e consistent with a virialized BLR 
(|Ulrich fc Horneiri996[) . 


1.3. The AGN STORM Project 

Given the importance of the UV emission lines in the 
photoionization equilibrium of the BLR gas and the prob¬ 
able differences between the geometry and kinematics of 
the high and low-ionization gas in the BLR, we have un¬ 
dertaken a large RM program in the UV using th e Gosmic 
Origins Spectrograph (GOS; iGreen et al.ll2012l) on HST 
{HST Program GO-13330), the AGN Space Telescope 
and Optical Reverberation Mapping (AGN STORM) 
Project, in the first half of 2014. The program was de¬ 
signed with certain specific goals in mind: 

1. Determine the structure and kinematics of the 
high-ionization BLR through observations of the 
variations in the G iv A1549, Lya A1215, N V A1240, 
SiivA1400, and HeiiA1640 emission lines 

2. Carry out simultaneous ground-based observations 
of (a) the high-ionization optical line HeiiA4686 
for direct comparison with HeiiA1640 and (b) 
the Balmer lines, particularly H/3A4861, to de¬ 
termine the structure and kinematics of the low- 
ionization BLR. Although the optical spectrum 
is extremely well-studied (iPeterson et al.l 20021 : 
iBentz et al.ll2007t l2010al IDennev et al.l 12010 , and 
references therein), simultaneous observations are 
necessary, as the dynamical timescale for the BLR 
in NGC 5548 is only a few years. 

3. Compare in detail the continuum variations in the 
UV (at ^ 1350 A) with those at other wavelengths 
fsee lEdelson et aHl2015l hereafter Paper 11) and in¬ 
fer the structure of the continuum-emitting region. 

The motivation for the UV/optical continuum compari¬ 
son is multifold: 


We note that three of these lines are actually doublets: 
NvAA1239,1243, Siiv A1394,1403, and C iv AA1548,1551. More¬ 
over, the He II feature is blended with O III] AA1661, 1665 and Si iv is 
blended with the quintuplet Oiv] AA1397.2, 1399.8, 1401.2, 1404.8, 
1407.4, where the second, third, and fifth transitions dominate. 
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1. Delays between continuum variations at longer 
versus shorter wavelengths have been de¬ 
tected or hinted at in a number of sources 


(e.g.. 

Wanders et al. 

119971: ICollier et alJ 

1998, 

2 OO 1 I: 

Peterson et al. 

119981: ISergeev et al.l 

2005; 

Cackett. Horne. & Winkleil 120071 IM“Hardv et al. 

2014 iShappee et al.l |2014jl. Such delavf 

can 


provide insight into the structure, geometry, and 
physics of the continuum-emitting region. 


2. Velocity-delay maps recovered using the UV con¬ 
tinuum as the driving light curve (Equation [IJ are 
expected to be of higher fidelity than those ob¬ 
tained from the optical continuum because the ob¬ 
servable UV is closer in wavelength to the ioniz¬ 
ing continuum (A < 912 A) that powers the emis¬ 
sion lines. The optical continuum is not only a 
slightly time-delayed version of the UV contin- 
uum, but it seems smoothed somewhat as well 
(|ShaDDee et al.l 120141: iPeterson et al.11201411 , which 
might make it difficult to recover detailed structure 
in the velocity-delay maps. 

Our HST program afforded a valuable opportunity for 
exploring AGN behavior at high time resolution for an 
extended period at wavelengths beyond those covered by 
our HST COS spectra. The HST program is the an¬ 
chor of a much broader AGN STORM project to address 
broader issues through observations across the electro¬ 
magnetic spectrum. This paper serves as the first in a 
series. 

Of special interest is the possibility of using short- 
timescale lags between variations in different continuum 
bands to map the temperatu re structure of the accre¬ 
tion disk. The Swift satellite (jGehrels et al.l[2n?)3 l is es¬ 
pecially suitable for such a study because of its broad 
wavelength coverage (hard X-ray through U-band) and 
ability to execute high-cadence observations over an ex¬ 
tended period of time. In Paper II, we present the re¬ 
sults of a four-month program of high-cadence (approxi¬ 
mately twice per day) multiwavelength observations with 
Swift. Additional papers in this series will describe high- 
cadence ground-based photometry from the near UV 
through the NIR. We will also present results from a pro¬ 
gram of ground-based spectroscopy that is similar in ca¬ 
dence to the HST GOS observations, but covers a some¬ 
what longer temporal baseline. Other additional papers 
will present results on the variable absorption features 
and on our efforts to decipher the broad emission-line 
variations and determine the structure and geometry of 
the BLR. 

In Section 2, we describe the observations and data 
processing, including a discussion of the program design 
and a complete description of how the standard data re¬ 
duction pipeline was modified to meet our stringent cali¬ 
bration requirements. We describe our initial data anal¬ 
ysis and results in Section 3, and in Section 4, we briefly 
discuss the first results from our program and place these 
results in the context of previous monitoring campaigns 
on NGC 5548. When necessary, we assume a ACDM cos¬ 
mology wit h Hg = 70km~| s~^ M pc~^, Dm = 0.28, and 
Da = 0.72 (|Komatsu et al.l[2011h . 

2. OBSERVATIONS AND DATA REDUCTION 


2 . 1 . Program Design 

RM is a resource-intensive activity that requires ob¬ 
taining high signal-to-noise ratio (S/N) homogeneous 
spectra at sufficiently high spectral resolution to resolve 
the gross kinematics of the BLR. Spectra must be ob¬ 
tained at a high cadence over a temporal baseline that is 
longer than the typical variability timescale of the AGN. 
Given the inherent risks of RM programs due to the un¬ 
predictability of AGN variability, it is essential that our 
experimental design assures a successful outcome, yet is 
as economical with observing time as possible. The first 
consideration is that each epoch of observation should 
require no more than one HST orbit per “visit” which 
restricts the integration time per visit to ~ 45 — 50 min¬ 
utes. This consideration limits us to relatively bright 
nearby Seyfert 1 galaxies. COS is clearly the instrument 
of choice for such a project, as it is a very sensitive, high 
spectral resolution spectrometer. Its native resolution 
(i? > 20000) is high enough to allow us to trade off res¬ 
olution and S/N in the data processing phase. In order 
to schedule the observatory efficiently, a cadence of one 
visit per day or longer is required. 

We therefore want to target an AGN that has a C iv- 
emitting region several light days in extent, and this 
requi res a source with logLx(1350A)/(ergss~^) A 43.5 
(e.g., iKaspi et al.l l2007ll . This led us immediately to 
select as a target the well-studied Seyfert 1 galaxy 
NGC 5548 (z = 0.017175). NGC 5548 is probably 
the best-studied AGN by RM, with historical optical 
spectroscopy extend ing as far back as the early 1970s 
(|Sergeev et al.l 1200711 . Importantly, it has never been 
known to go into a “dor mant state,” as obse rved recently 
in the case of Mrk 590 (jPennev et al.l [201^ . that would 
preclude a successful reverberation campaign and, histor¬ 
ically, self -absorption in the UV resona nce lines has been 
minimal (ICrenshaw fc Kraemeii Il999f). although strong 
absorption appeared in 2013 (|Kaastra et al.ll2014ll . 

The remaining adjustable parameter is the dura¬ 
tion of the campaign. We investigated this using 
Monte Carlo simu lations similar to those described by 
iHorne et al.l (I2004D . Using recent dev elopments in sta¬ 
tistically modeling AGN light curves (iKellv et al.l 1^0091 
iKozlowski et al.l 120101 IMacLeod et al.l 120100 , we can 
make very robust models of the expected continuum 
behavior of NGC 5548. Quasar light curves are well- 
described by a stochastic process, the damped ran¬ 
dom walk. The process is described by an amplitude 
tr and a damping timescale rd, which for NGC 5548 
in the optical are measured to be cr = 0.891 q 2 o ^ 

10 “^^ ergs s“^cm“^ A ^ and Td = 77+34 days, respec¬ 
tively (IZu. Kochanek. fc Petersonll2011ll . We used these 
measured properties of NGC 5548 to simulate the con¬ 
tinuum variations; this is a conservative choice as the 
UV continuum can be expected to show both higher am¬ 
plitude and shorter time-scale variations, both of which 
are an advantage. We then convolved the artificial 
light curves with model velocity-delay maps for several 
lines to prov i de an artificial spectrum. As described by 
iHorne et al.l (|2004D . we adopted a BLR model with an 
extremely challenging velocity-delay map for these sim¬ 
ulations, a Keplerian disk with a single two-armed spiral 
density wave. While this is unlikely to be the actual 
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AGN BLR geometry, it provides a challenging test: if we 
can recover such a complex velocity-delay map correctly, 
then we can certainly hope to recover others of com¬ 
parable complexity and would have no difficulties with 
geometries lik e those that have been recovered for opti¬ 
cal lines le g.. iBentz et ani2ninbl: iPancoast et al.l[20121 : 


iGrier et al.ll2013l: iPancoast et alJl2014 ). 

We modeled the emissivity and response of each line 
realistic ally using a grid o f photoionization equilibrium 
models (|Horne et al.ll20d^ . We sampled the artificial 
spectra to match our proposed observations, including 
noise. We then modeled the artificial spec tra to recover 
the velocity-delay maps using MEMECHO (|Horn i IT^ 
iHorne et al.ll2004r) . Simulations based on characteristics 
of previous RM experiments yield velocity-delay maps 
with noise levels similar to those obtained from the actual 
dat a, demonstrating th e verisimilitude of our simulations 
(see IHorne et al.l[2004 for examples). 

The goal of our simulations was to determine the min¬ 
imum duration program that would allow us to recover a 
velocity-delay map with a high probability of success. 
For COS-like observations (in terms of S/N per visit 
and spectral resolution), our initial simulations indicated 
that reliable velocity-delay map recovery for a strong line 
(e.g., C iv) required between 130 to 200 days. A finer grid 
of models showed that in 10 of 10 simulations, a high- 
fidelity velocity-delay map was recovered after 180 days, 
which was thus adopted as the program goal. A much 
longer program at this sampling rate would in any case 
be precluded by the accessibility of the target to HST. 

Because of the long duration of the proposed program, 
we also considered the possible impact of losses of data 
due to instrument or spacecraft safing events. Short sat¬ 
ing events occur frequently enough that we needed to as¬ 
sess their impact. Based on the record for HST and GOS 
in Gycles 17-20, there might be two spacecraft events 
that lose 2-3 days each and one GOS event that loses 2 
days over a stretch of 180 consecutive days. By repeat¬ 
ing a subset of our simulations, we found that losses of 
such small numbers of observations would have no im¬ 
pact on our ability to recover the velocity-delay maps. 
The simulations also allowed us to assess the impact of 
early termination of our experiment due to a major fail¬ 
ure. If a program was terminated at ^ 100 days, the 
probability that the data would yield a useful (but not a 
detailed) velocity-delay map would be ^ 50%. However, 
a program as short as 75 days would have a very low 
probability (^ 10%) of success. 

The key to a successful RM campaign is that it must be 
long enough that favorable continuum variability char¬ 
acteristics become highly probable. That this is essen¬ 
tially guaranteed to happen during a 180-day experiment 
played a major role in selecting NGC 5548 as our target. 

2 .2. COS Observations 


Observations were made in single-orbit HST GOS visits 
approximately daily from 2014 February 1 through July 
27. Of the 179 scheduled visits, 171 observations were 
executed successfully and 8 were lost to safing events 
or target acquisition failures (very close to the expected 
number of losses). 

In each visit, we used the G130M and G160M gratings 
to observe the UV spectrum over the range 1153-1796A 
in four separate exposures. Exposure times were selected 


to provide S/N ^ 100 when measured over velocity bins 
of ~ 500kms~ . During each visit, we obtained two 
200-second exposures with G130M centered at 1291A 
and 1327 A and two 590-second exposures with G160M 
centered at 1600 A and 1623 A. 

The GOS far-ultraviolet detector is a windowless, 
crossed delay-line microchannel plate stack that is sub¬ 
ject to long-term charge depletion. To extend the use¬ 
ful lifetime of the detector, we positioned the spectrum 
so that bright geocoronal airglow lines (e.g., LyaA1215) 
and AGN emission lines (e.g., redshifted Lya) would not 
always fall on the same area of the detector. First, we al¬ 
ternated the target acquisition between the G130M/1291 
and the G130M/1327 configurations. The G130M/1327 
configuration is then followed by a G130M/1327/FP- 
POS=3 exDOSurj^. and by a G130M/1291 exposure al¬ 
ternating among FP-POS=l, 2, and 4. The G130M/1291 
configuration is instead followed by a G130M/1291/FP- 
POS=3 exposure, and by a G130M/1327 exposure al¬ 
ternating among FP-POS=l, 2, and 4. Second, we al¬ 
ternated the FP-POS for the G160M/1623 exposure be¬ 
tween FP-POS=l and FP-POS=2. We could not vary 
the settings for the G160M/1600 and G160M/1623 fur¬ 
ther because we needed to ensure the coverage of the 
entire wavelength range while keeping the detector gap 
from falling on the redshifted Giv A1549 emission line. 

Finally, we used four additional orbits to improve our 
understanding of the GOS flux calibrations (see Section 
Ell). During these additional visits, we observed two of 
the standard stars (WD 0308-565 and W D 1057-1-719) 
empl oyed to obtain sensitivity functions (|Massa et al.l 
1201411 at the same detector locations we used for the re¬ 
verberation program. The observations were taken using 
all the instrument configurations employed in our pri¬ 
mary observing program. 

2.3. Data Reduction 

We used the CalCDS pipeline v2.21 for the bulk of our 
data processing. The absolute flux calibration of the 
GOS reduction pipeline is reported to be accurate to 
5% and the relative flux calibrat ion is good to bet¬ 
ter than ^ 2% (|Holland et ahllMIll . We are primarily 
interested in the quality of the relative flux calibration 
as we are looking for very small-scale variations on short 
timescales; we need the fluxes to be stable and repeat- 
able across the spectrum. We found, however, that there 
were local variations in the precision of the fluxes that 
necessitated improvements. 

To produce a final dataset with a flux calibration that 
is everywhere precise at the 2% level, we refined the exist¬ 
ing calibration reference files and applied a post-CalCOS 
pipeline to further process the data. The main areas of 
improvement include refinements to the dispersion solu¬ 
tion, fixed-pattern noise mitigation, the sensitivity func¬ 
tion, and the time dependent sensitivity (TDS) functions, 
as outlined below. The final data product consists of one 
combined spectrum per grating per day. Airglow emis¬ 
sion lines (OiAA1302.2,1306 and NiAA1199.5,1200.7) 
were filtered from the data by removing events detected 

43 FP-POS values refer to small displacements of the spectrum 
on the detector in the dispersion direction in order to minimize the 
effects of fixed-pattern noise. 
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when HST was in daylight. The spectra were further 
binned by 4 pixels in order to increase the S/N per spec¬ 
tral element of the AGN continuum. This binning still 
results in two binned pixels per COS resolution element. 

2 . 3 . 1 . Dispersion solution 

The COS wav elength solution has a quoted uncertainty 
of ^ 15kms“^ (jHolland et al.l 1^01411 . As the STIS un¬ 
certainty is < 5kms“^ (jHernandez et al.1120141 1. we re¬ 
fined the dispersion solutions for our COS dataset using 
previous observations of NGC 5548 taken with the STIS 
E140M/1425 mode in 1998 (PID 7572, PI: Kraemer). To 
accomplish this, we cross-correlated the line profiles of 
strong interstellar medium absorption features between 
each COS observation and the STIS reference spectrum. 
We used 19 interstellar absorption features, ranging from 
SiiiA1190 at the short-wavelength end to AliiA1670 at 
the long-wavelength end. A linear correction to the ini¬ 
tial wavelength solution was then computed across each 
detector segment and applied directly to the extracted 
spectra. With this correction, measurements of the root 
mean square (rms) of the residual offsets decreased from 
~ 15kms“^ to < 6kms“^. 

2 . 3 . 2 . Fixed-pattern noise 

The standard reference files used in the CalCOS pipeline 
correct for only the most prominent fixed-pattern noise 
features such as the quantum-efficiency gridwires, low- 
order res ponse variations , and large geometric distortion 
artifacts (jElv et al.ll20lH) . Usually, users combine mul¬ 
tiple FP-POS positions to smooth over the remaining 
features. However, this was not possible for our dataset, 
since only a single FP-POS setting was used for each cen¬ 
tral wavelength setting in each orbit (see Section 12.2|) . 

To correct these features to a higher degree we 
derived one-dimensional pixel-to-pixel flats (“p-flats”). 
These flats were produced by combining normalized, 
high signal-to-noise ratio white dwarf spectra in detec¬ 
tor sp ace, following the method described by lElv et all 
(1201111 . The white dwarf spectra used were taken as part 
of the HST/CAL program 12806 (PI: Maasa) and used 
the same detector locations as used for the NGC 5548 
datasets. 

To test the effects of our p-flat correction, we com¬ 
bined the 171 spectra reduced both with and without 
the application of the p-flats. The S/N per pixel in 5 A 
continuum regions increased from ^ 75 to ^ 80 for the 
G130M grating, and from ^ 60 to ^100 for the G160M 
grating through the removal of small localized flux cali¬ 
bration errors by the p-flat correction. The improvement 
for the G130M combined spectrum is less dramatic be¬ 
cause (a) we rotated among the four FP-POS settings, 
and (b) the G130M grating disperses more widely in the 
cross-dispersion direction, and thus intrinsically averages 
the fixed-pattern noise over a larger area of the detector. 

2 . 3 . 3 . Sensitivity functions and TDS 

The COS flux calibration is done in two steps: (a) 
derivation of static sensitivity functions and (b) charac¬ 
terization of the tim e evolution of the sen sitivity through 
the TDS correction ((Holland et al.l[^14H . Thanks to the 
existing calibration program that monitors the TDS vari¬ 
ations (PID 13520), we had bi-monthly observations of 


the standard star WD 0308-565 for 3 out of the 4 cen¬ 
tral wavelength settings we are using in our program 
(G130M/1291, G130M/1327 and G160M/1623). Stan¬ 
dard star data were obtained in 2014 February, April, 
June, and August. By analyzing these calibration data, 
together with the data c ollec ted during our additional 
calibration orbits fSection [^1^ . we verified that both the 
static and time-dependent response functions vary more 
with instrument configuration than currently modeled by 
the CalCOS pipeline. While CalCOS assumes that both 
the sensitivity function and the TDS correction vary only 
as a function of wavelength, we were able to improve the 
relative flux calibrations and reach our required level of 
precision by (a) obtaining sensitivity functions individu¬ 
ally for each configuration (one function per wavelength 
setting per FP-POS per detector segment), and (b) com¬ 
puting the TDS correction individually for each wave¬ 
length setting observed as part of the routine calibration 
program. 

We estimate the quality of the flux calibration by in¬ 
spection of the fractional residuals /res of the calibrated 
standard star spectra and their respective CalSPEC stel¬ 
lar model (the same models employed by the standard 
pipeline reference file), 

r _ /wD ,obs — /Model 

/res — 7 ) (^) 

/Model 

where both /wD.obs and /Model are binned over 1A using 
a boxcar filter in order to increase the S/N per spectral 
element. While the visual inspection of the residuals as a 
function of wavelength allows us to identify and correct 
for local biases, we use the mean value of the distribu¬ 
tion of the residuals as an indicator of a global bias in 
the calibration. The flux calibration uncertainty is an 
estimate of the limit of the stability of the flux calibra¬ 
tion at a given time. However, since the overall instru¬ 
ment sensitivity evolves with time, and our final spectra 
are obtained from the combination of multiple settings 
for each grating, we conservatively define the fractional 
precision error Jp for each grating as the maximum frac¬ 
tional uncertainty computed for any of the wavelength 
settings. 

The new sensitivity functions were derived from spec¬ 
tra of the standard star WD 0308-565 for the G130M 
settings, and of WD 1057-1-719 for the G160M settings 
(PID 12806). While one individual sensitivity function 
per detector segment characterizes the full grating (data 
from different settings are averaged together) in CalCOS, 
we built one independent sensitivity function for each 
wavelength setting and FP-POS used in our program. 

By comparing the bi-monthly WD 0308-565 data, we 
found that residuals with respect to the stellar models 
were greatly reduced if the TDS corrections were com¬ 
puted individually for each of the wavelength settings 
(G130M/1291, G130M/1327 and G160M/1623), instead 
of averaging the data over multiple modes. Additional 
improvements were obtained by increasing the number 
of time intervals over which the TDS trends are com¬ 
puted and by redefining the wavelength ranges used in 
the analysis. Unfortunately, there are insufficient cali¬ 
bration data for the longest wavelengths in the G160M 
spectra, so we were forced to truncate these spectra at 
1750 A. 
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In spite of these improvements, the available data did 
not allow us to conduct any tests on the remaining set¬ 
ting (G160M/1600). This setting is particularly impor¬ 
tant for our scientific goals since it includes most of the 
broad Civ emission line fsection ll.3p . Moreover, since 
all the calibration data for TDS monitoring purposes 
are obtained only using FP-POS=3, they did not al¬ 
low us to test for any residual dependence of the TDS 
correction on FP-POS configuration. These are the two 
reasons that motivated us to request further calibration 
data (see Section [52]) • These data for WD 0308-565 and 
WD 1057-1-719 collected in 2014 September allowed us 
to derive an independent set of sensitivity functions. By 
comparing the new sensitivity functions with the orig¬ 
inals, appropriately corrected for time evolution of the 
TDS, we were able: 

1. To identify the best possible TDS correction attain¬ 
able for the G160M/1600 setting with the current 
TDS calibration data. The current CalCOS TDS 
correction for this configuration was obtained from 
combining both G160M/1577 and G160M/1623 
data. Although this correction is not ideal, it min¬ 
imizes both the global bias and the flux calibration 
uncertainty when compared to the TDS corrections 
obtained individually from either the G160M/1577 
or the G160M/1623 settings. 

2. To verify that the TDS correction does not vary 
strongly with FP-POS settings. While compar¬ 
ing the residuals with respect to the stellar models 
shows structure unique to each FP-POS setting, 
the level of the local biases is such that both the 
global bias and the flux calibration uncertainty can 
be considered stable for each wavelength setting 
(the maximum deviation in the width of the resid¬ 
ual distribution is ~ 0.2%). 

With the new flux calibration and TDS characteriza¬ 
tion, our global biases are consistent with zero and the 
overall fractional precision is dp ~ 1.1% and ~ 1.4%, re¬ 
spectively, for the G130M and G160M settings, compared 
to dp ^ 1.4% and ~ 3.5% for the standard pipeline. We 
intend to make these improvements available to other 
GOS users. 

2.3.4. Sensitivity offsets for the final week of data 

During the final week of the observing campaign, the 
operating high voltage (HV) for one of the GOS detec¬ 
tor segments was i ncreased to combat the negative ef¬ 
fects of “gain-sag” ([Sahnow et al.ll2MTI) . While neces¬ 
sary to provide well-calibrated data, this HV change also 
has the effect of introducing small changes in the de¬ 
tector response. Using WD 0308-565 observations from 
our calibration orbits taken at the same HV as the rest 
of the campaign and contemporaneous TDS monitoring 
observations of the same calibration target (taken at the 
increased HV), we were able to estimate a HV bias cor¬ 
rection for the G130M grating from a direct comparison 
of the spectra. This bias was measured to be 1% with¬ 
out any detectable dependence on wavelength or cen- 
wave setting. Unfortunately, the same procedure could 
not be done for the G160M grating as our additional or¬ 
bits and the TDS observations used different standard 
stars (WD 0308-565 instead of WD 1057-1-719). Instead, 


the bias estimate for this grating was obtained by an¬ 
alyzing the time evolution of the mean flux in overlap¬ 
ping regions of G130M/1327 at the lower HV setting and 
G160M/1600 at the higher HV setting. This analysis 
gave a plausible estimate of a 1% bias. However, with 
such a limited amount of data at the lower HV and a 
narrow overlapping wavelength range, the estimate lacks 
the accuracy of the G130M bias estimate. 

3. DATA ANALYSIS 
3.1. Mean and RMS Spectra 

For an initial look at the spectral variations, we define 
G130M and G160M mean spectra as 

_ 1 ^ 

= (3) 


where Fi is the fth spectrum of the series of V = 171 
spectra. Similarly, the rms residual spectrum (hereafter 
referred to simply as the RMS spectrum) is defined as 


S{X) = 


N -1^ 

i=l 


1/2 


[F{x)-F{x)y 


( 4 ) 


The RMS spectrum is especially useful as it isolates the 
variable part of the spectrum; constant components dis¬ 
appear, though sometimes small residuals are visible in 
the case of strong features. 

The statistical uncertainty in the mean spectra is 


W — AT 



1/2 


^Fi (^) 


( 5 ) 


where aPi is the error spectrum of the fth spectrum in 
the series. 

The total uncertainty in the mean spectra consists of 
this statistical uncertainty and our estimate of the frac¬ 
tional uncertainty in precision as described above, which 
amounts to dp(G130M) Ri 1.1% and (5p(G160M) « 1.4%. 
To determine the total uncertainty, the statistical un¬ 
certainty (Equation |S]) and the uncertainty in precision 
(dpF(A)) are added in quadrature. The mean and RMS 
spectra for the G130M and G160M settings are shown in 
Figures [T] and m respectively. 

The RMS spectrum resulting from Equation (|3|) com¬ 
bines both the intrinsic varia bility and the var iance 
due to nois e , as discussed by iPark et ^ (|2012t) and 
iBarth et al.l (j2015l l. In order to isolate the RMS spec¬ 
trum of the intrinsic variations tro; we model the distri¬ 
bution of the residuals of each pixel about the mean. The 
combined statistical and systematic noise in each pixel 
of spectrum i is thus {af + (ApFi)^)^/^. Assuming that 
the flux measurement errors and the intrinsic variations 
arise from independent Gaussian random processes, we 
find maximum likelihood estimates for the optimal aver¬ 
ages, F^, and (Tq by minimizing 


N 


-2 InL (F;,,cro) = X^+Eln + {5pF t) , (6) 


The unweighted average F in Equation |3] is formally distinct 
from the optimal average E/j, though practically they are indistin¬ 
guishable for these data. 
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Figure 1. Spectra obtained with the G130M grating. The shaded areas show the integration regions defined in Table HI . Top panel: 
the mean spectrum, as defined in Equation is shown as a black solid line. The spectral model described in Section IT5I is shown in 
red. The deep trough centered at 1215 A is Galactic Lya absorption and the strong narrow emission line at the center of this trough is 
geocoronal Lyo emission. The narrow Galactic absorption lines, although generally saturated, are never black at line center because the 
thermal broadening of this cold gas is still far below the resolution of COS. Botto m pa nel: the RMS spectrum, as defined in Equation 
is shown as a black solid line, while the intrinsic RMS spectrum (jq (see Section [311 is shown in blue. In grey we show the total error 
on the mean, which is the statistical error (Equation O combined in quadrature with the fractional error in precision, = 1.1% for the 
G130M spectra. Note the difference in the flux scale between the two panels. 


where 


X 




(7) 


and SpFi is the precision of the ith spectrum. This es¬ 
timate of the intrinsic RMS spectrum (cro) is also shown 
in Figures [T] and m 


3.2. Integrated Light Curves 

The next step in our initial analysis is to produce light 
curves for the continuum and emission lines. At this 
stage, our goal is to make simple measurements from the 
reduced spectra, introducing as few assumptions as pos¬ 
sible. All flux measurements are performed on spectra in 
the observed frame. We have not corrected the spectra 
for Galactic extinction in order to facilitate the clean¬ 
est comparison with other measurements to be reported 
elsewhere in this series of papers (e.g., broad-band pho¬ 
tometry). 

There are bad pixels throughout the spectrum, and 
their location and severity change with time, instrument 
settings, and airglow subtraction (e.g., if a spectrum is 
taken entirely in orbital bright time, the flux in the air- 


glow windows is set to zero and the pixels are flagged 
as bad pixels). To prevent the introduction of artificial 
variations in the relative flux estimates, bad pixels are 
masked throughout the dataset. This means that if a 
pixel is bad in any of the visits, the pixel is masked out 
in each of the 171 spectra. We further mask Galactic 
Lya absorption and airglow region. Integration ranges 
(listed in Table [1} were chosen using the mean spectra in 
Figures [Hand [2] as a guide. Gontinuum ranges are chosen 
to be as uncontaminated as possible by absorption lines 
and broad emission-line wings. In the case of overlap¬ 
ping emission lines (e.g.. Civ and Hen), the boundary 
wavelength corresponds to the wavelength at which the 
fluxes of the two lines are comparable. We do not mask 
absorption lines at this stage in our analysis. We are un¬ 
able to cleanly separate N v and Lya using this simple 
procedure. 

Continuum fluxes are measured as the weighted mean 
of the flux density in the integration region, with weights 
equal to the inverse of the variance. 
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Figure 2. Spectra obtained with the G160M grating. The shaded areas show the integration regions defined in Tabl e fT|l . Top panel: the 
mean spectrum, as defined in Equation 0, is shown as a black solid line. The spectral model described in Section IsTsi is shown in red. 
The narrow Galactic absorption lines, although generally saturated, are never black at line center because the thermal broadening of this 
cold gas is still far below the resolution of COS. Botto m pa nel: the RMS spectrum, as defined in Equation is shown as a black solid 
line, while the intrinsic RMS spectrum <jq (see Section 13.111 is shown in blue. In grey we show the total error on the mean, which is the 
statistical error (Equation O combined in quadrature with the fractional error in precision, 5p = 1.4% for the G160M spectra. Note the 
difference in the flux scale between the two panels. 


Table 1 

Integration Limits for Light Curves 


Emission 

Component 

Integration 

Limits 

Shortward 
Continuum Region 

Longward 
Continuum Region 

Fa (1367A) 
Lya A1215 
SiivAllOO*" 
CivA1549 
Hell A1640'" 

1364.5-1369.5 

1201.0-1255.0 

1405.0-1455.0 

1520.0-1646.0 

1647.0-1740.0 

1155.0-1160.0 

1364.5-1369.5 

1475.0-1482.0 

1475.0-1482.0 

1364.5-1369.5 

1460.0-1463.5 

1743.0-1749.0 

1743.0-1749.0 

Note. — All 

regions are in 

the observed frame (A) 



^ Integration range also includes Oiv] A1402. 
^ Integration range also includes O ill] A1663. 


where Wi = as in Equation ([5]). Statistical uncer¬ 
tainties c omputed by Ca lCOS are corrected for low counts 
following I Gehrel j (|1986ll . Statistical uncertainties on the 
mean fluxes are obtained through standard error propa¬ 
gation, 





- 1/2 


(9) 


In all cases, bad pixels are excluded from the computa¬ 


tion. 

Emission-line fluxes are measured as the numerical in¬ 
tegral of the emission flux above a locally defined contin¬ 
uum defined by the relatively featureless windows given 
in Table [TJ To estimate the local continuum underneath 
the line we performed a linear fit of the continuum 
flux in the selected regions. The linear local continuum 
is then subtracted from the emission component, again 
masking bad pixels. The line flux is numerically inte¬ 
grated over the integration limits given in Table [T] us- 
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ing Simpson’s method. We do not interpolate over bad 
pixels. We note, however, that the difference between 
integrating over the bad pixels and computing the in¬ 
tegral excluding them is < 0.1%. Statistical errors are 
computed numerically by creating A’sampie = 5000 real¬ 
izations of the line flux and the underlying linear contin¬ 
uum. The flux Fx is randomly generated from a Gaussian 
distribution having mean equal to the flux of the spectral 
element and width a equal to the statistical error on the 
flux. For the linear continuum, we generate A’sampie fits 
having a mean equal to the best fit values and covari¬ 
ance equal to their covariance matrix. For each realiza¬ 
tion, a line-flux estimate is then obtained by subtracting 
the linear continuum and by performing the numerical 
integration of the residuals. Confidence levels (lu) are 
finally obtained from the distribution of the A^sampie line 
fluxes. When the error bars are asymmetric, we adopt 
the larger error as the statistical error associated with 
the integrated flux. 

As noted above, we adopt as the fractional error in 
precision Jp = 1.1% and = 1.4% for the G130M and 
G160M settings, respectively. This is added in quadra¬ 
ture to statistical error of the integrated fluxes. The error 
in the precision dominates throughout the G160M spec¬ 
tra and in the G130M spectra as well, except at wave¬ 
lengths shortward of ~ 1180 A longward of ~ 1425 A, 
and in the core of the Lya complex. 

The final continuum and emission-line light curves are 
listed in Tableland shown in Figure |31 The light curve 
statistics are given in Table |3l The average interval be¬ 
tween two consecutive observations is (At) = 1.0 days 
with an rms cr(At) = 0.3 days. The median interval be¬ 
tween observations is Aimed = 1-0 days. The largest gaps 
between consecutive observations are three days (on two 
occasions) and two days (on six occasions). 

3.3. Time-Series Measurements 

Certain simplifying assumptions underlie the RM tech¬ 
nique. Most time-series analyses start with the assump¬ 
tion that the emission-line light curves are simply scaled, 
time-delayed, and possibly smoothed versions of the con¬ 
tinuum light curve. Inspection of the light curves in Fig¬ 
ure [3] suggests that this is an entirely reasonable assump¬ 
tion for the first half of the campaign. However, approx¬ 
imately halfway through the campaign, the emission-line 
response becomes more complicated. Between approxi¬ 
mately HJD2456780 and 2456815, the emission-line light 
curves are either fiat (Lya, Hell) or decreasing (Civ) 
while the continuum is slowly rising. Moreover, the in¬ 
tensity ratio between the last two strong peaks in the con¬ 
tinuum light curve at around HJD2456820 and 2456840 
seems to be almost inverted in the lines, with the sec¬ 
ond peak being stronger than the first one (especially in 
Civ). There is also a small event in the continuum light 
curve around HJD2456785 that does not appear to have 
counterparts in the emission-line light curves. The line 
light curve that seems to best trace the continuum is the 
Hell light curve, which is sensitive to the continuum at 
energies above 4 Ryd. It is also the only strong line in 
the COS spectra that is neither a resonance line nor self- 
absorbed. Moreover, it is the line that arises closest to 
the continuum source, as we will show below. 

Because of the changing character of the emission-line 
response, for our initial analysis we measure emission¬ 


line lags (a) for the entire data set and (b) for subsets 
that divide the data into two separate halves of 85 obser¬ 
vations each. The first subset, which we will refer to as 
“TI,” runs from HJD2456690 to 2456780 and the second 
subset, “T2,” runs from HJD245678I to 2456865. 

We first measured the emission-line lags relative to 
the continuum variations by cross-correlation of the 
light curves. We used the interpolation c ross-correlation 
ICCF) method as implemented by iPeterson et al.l 
20041) . In this method, uncertainties are estimated us¬ 
ing a mo del-independent Monte Carlo method referred 
to as “flux randomization and random subset selection 
(FR/RSS).” For each realization, N data points are se¬ 
lected from a light curve with N independent values, 
without regard to whether or not any particular point 
has been previously selected. For data points selected n 
times in a given realization, the flux error associated with 
that data point is reduced by a factor of The flux 

measured at each data point is then altered by adding 
or subtracting a random Gaussian deviate scaled by the 
flux uncertainty ascribed to that point. Each realiza¬ 
tion yields a cross-correlation function that has a max¬ 
imum linear correlation coefficient r„iax that occurs at 
a lag Tpeak- We also compute the centroid Tcent of the 
cross-correlation function using all the points near Tpeak 
with r(T) > O.Srmax- Typically a few thousand realiza¬ 
tions are used to construct distribution functions for the 
ICCF centroid and peak. We adopt the median values of 
the cross-correlation centroid distribution and the cross¬ 
correlation peak distribution as our lag measurements. 
The uncertainties, which are not necessarily symmetric, 
correspond to a 68% confidence level. In general, Tcent is 
found to be a more reliable indicator of the BLR size than 
T'peak, though we record both. The ICCF measurements 
of Tpeak and Tcent for the four strongest UV emission lines 
are given in the second and third columns, respectively, 
in Table m 

We have also estimated emission-line lags using 
JAVELIN, which is an improve d version of SPEAR 
(|Zu, Kochanek, fc Peterson! I2011D . JAVELIN assumes 
that the emission-line light curves are shifted and 
smoothed versions of the continuum light curve (as with 
the ICCF analysis), w here the continuum is modeled as a 
damped random walk (iKellv et alll2009l iKozlowski et al.l 
l20I0t iMacLeod et al.ll20I0D with uncertainties deter¬ 
mined using the Markov Chain Monte Carlo method. We 
model the full dataset, and each line light curve was run 
independently with the continuum. The results, given in 
column (4) of Table |4j are in good agreement with the 
ICCF analysis, as expected. 

In columns (5) and (6) of Table SI we also give the 
ICCF centroid values for the TI and T2 subsets. We 
also show the ICCFs for the entire sample and the TI 
and T2 subsamples in Figure SI In general, the lags for 
the TI subsample have the smallest uncertainties and the 
ICCFs have the largest peak correlation coefficients Tmax, 
as expected. The T2 subset, on the other hand, yields 
lags with larger uncertainties and ICCFs with lower val¬ 
ues of r„iax (indeed, much lower in the case of Siiv and 
Civ), again as expected from visual inspection of the 
light curves. The T2 lags are also larger than those from 
TI, probably only in small part because the continuum 
is on average brighter (by ^ 15% on average) during the 
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Table 2 

Continuum and Emission-Line Light Curves 



G130M 




G160M 


HJD^ 

Fx (1367A)® 

F(Lya)'^ 

F(Si !¥)== 

HJD'^ 

F{C !¥)<= 

F(Heii)‘= 

6690.6120 

34.27 ± 0.64 

39.66 ± 

0.47 

4.04 ±0.26 

6690.6479 

53.24 ±0.79 

6.92 ±0.31 

6691.5416 

35.45 ±0.65 

39.88 ± 

0.48 

4.47 ±0.30 

6691.5760 

53.06 ±0.79 

6.99 ±0.34 

6692.3940 

37.71 ± 0.67 

39.88 ± 

0.48 

4.83 ±0.27 

6692.4084 

53.30 ±0.80 

6.51 ±0.35 

6693.3237 

38.14 ± 0.68 

39.22 ± 

0.47 

4.19 ±0.28 

6693.3380 

53.08 ±0.80 

6.64 ±0.36 

6695.2701 

40.94 ± 0.71 

39.52 ± 

0.47 

3.92 ±0.29 

6695.3145 

53.09 ±0.81 

7.36 ±0.35 

6696.2459 

44.25 ±0.75 

39.49 ± 

0.48 

3.72 ±0.29 

6696.2602 

52.76 ±0.80 

7.25 ±0.36 

6697.3080 

45.30 ± 0.75 

40.16 ± 

0.49 

4.38 ±0.30 

6697.3223 

53.77 ±0.82 

8.00 ±0.36 

6698.3041 

48.27 ± 0.79 

40.04 ± 

0.48 

4.14 ±0.30 

6698.3184 

55.40 ±0.83 

8.75 ±0.36 

6699.2338 

45.80 ± 0.76 

41.43 ± 

0.51 

4.80 ± 0.35 

6699.2481 

55.65 ±0.84 

8.77 ±0.37 

6700.2299 

46.00 ± 0.76 

41.13 ± 

0.50 

4.37 ±0.30 

6700.2442 

55.13 ±0.83 

7.73 ±0.38 

6701.3588 

47.46 ±0.78 

41.75 ± 

0.50 

4.52 ±0.33 

6701.3731 

54.82 ±0.83 

8.41 ±0.38 

6702.1557 

47.74 ± 0.78 

41.98 ± 

0.51 

4.41 ± 0.34 

6702.1700 

55.64 ±0.84 

8.72 ±0.39 

6703.1518 

47.56 ±0.78 

42.33 ± 

0.51 

4.70 ± 0.32 

6703.1661 

55.63 ±0.84 

9.27 ±0.37 

6705.3432 

45.77 ± 0.76 

43.80 ± 

0.53 

5.57 ±0.33 

6705.3575 

57.85 ±0.86 

9.31 ±0.34 

Note. — Full table is given in the published version. 

Integrated light curves in 

the observed 

frame. Flux uncertainties include both statistical and systematic errors. 

® Midpoint of the observation (HJD — 2450000). 

^ Units of 10“^® erg s”^ cm“^ A“l. 

Units of 10“^® erg s“^ cm~^. 



Table 3 

Light Curve Statistics 

Emission 

Mean and 

Mean 

Maximum 

Minimum 


Component 

RMS Fiux 

Fractional Error 

TP ^ 

-f^var 

Fiux 

Fiux 

P b 

-O-max 

(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

Fx (1367A)'= 

42.64 ±8.60 

0.017 

0.201 

64.74 

21.87 

2.96 ±0.08 

F{Lya)d 

41.22 ±2.71 

0.012 

0.065 

46.84 

34.31 

1.37 ±0.02 

F(Siiv)d 

4.62 ±0.55 

0.065 

0.099 

6.00 

3.23 

1.86 ±0.15 

F(Civ)‘l 

53.28 ±3.91 

0.015 

0.072 

62.97 

47.23 

1.33 ±0.03 

F(Heii)‘i 

7.93 ± 1.13 

0.046 

0.135 

10.62 

5.47 

1.94 ±0.12 


Note. — Light curves statistics are in the observed frame. 
^ Excess variance, defined as 


(F) 


( 10 ) 


where cr is the RMS of the observed fluxes (column 2), S is t he mean statistical uncertaint y 
(column 3 times {-F}), and (F) is the mean flux in column (2) (|Rodrfauez-Pascual et ^1199711 . 
^ Ratio between maximum and minimum flux. 

^ Units of 10“^^ erg s“^ cm“^ A“^. 

^ Units of 10“^^ erg s“^ cm“^. 


Table 4 

Emission-Line Lags 


Emission Line 

^peak 

T”cent^ 

'TJAVELIN^ 

'^cent,Tl^ 

'^cent,T2^ 

Lyo 


6 19+0-29 

5-80t“ls 

5 90+0-20 
0.yU_Q 29 

IJ yq-|-0.76 
‘ • ''^-0.57 

Si IV 


5 44+0-70 

5-94lS:g 

4 99+0-76 

'T 99 -I-I .33 
' ‘^^-1.06 

Civ 

^•^- 0.6 

c oo-tO.44 
o.oo_Q 4 g 

4 59+0-68 
^•oy_0.42 

4 61+0-36 

^•^^-1.04 

Hen 

2 4+0-2 

2 c^n+0.34 
^•0'9 _o.31 

2 42+0-67 

9-11 -1-0.43 
^•-LJ-_o.38 

q 07 - 1 - 0.71 

■J-O'-0.58 


Note. — Delays measured in light days in the rest frame of NGC 5548. 
^ Complete dataset: 171 visits 
^ T1 dataset: visits 1-85 
^ T2 dataset; visits 86-171 
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second half of the campaign so the BLR gas that is most 
responsive to continuum changes is farther away from the 
central source. 

3.4. Velocity-Binned Results 

As noted earlier (Section 12.ip . this program was de¬ 
signed to recover kinematic information about the BLR 
by resolving the emission-line response as a function of 
radial velocity. This will be explored in detail in sub¬ 
sequent papers in this series. Here we carry out a sim¬ 
ple preliminary analysis intended to show only whether 
velocity-dependent information is present in the data. 
We isolate t he L yg and Civ profiles as described ear¬ 
lier ISection 13.21) but then integrate the fluxes in bins of 
width 500 km s“^, except in the shortward wing of Lya 
(—10,000 km s~^ < AH < —7000 km s~^) where we use 
1000 km s“^ bins on account of the low flux in the blue 
wing of this line. 

We show the ICCF centroids for each velocity bin in 
the C IV emission-line profile in the bottom panel of Fig¬ 
ure [S] For the two subsets as well as the entire dataset, 
we see that there is a clear ordered structure in the kine¬ 
matics. We cannot infer much from such a simple analy¬ 
sis, of course, because we cannot accurately characterize 
a complex velocity field with a single number. It is re¬ 
assuring, however, that the general pattern is similar to 
what has been seen in other objects and is qualitatively 
consistent with a virialized region (i.e., the high veloc¬ 
ity wings respond first). The middle panel of Figured 
shows the RMS spectra for the entire dataset and the 
two subsets. 

In the bottom panel of Figure [6l we show the velocity- 
binned ICCF centroids for the Lya emission line. Again, 
a clear pattern emerges, as the lags in each velocity bin 
are highly correlated with those of adjacent bins. How¬ 
ever, the pattern that emerges is unlike what is seen in 
CiV; the largest lags are at intermediate velocities and 
the lags decrease toward line center. However, given the 
severe blending and strong absorption, detailed model¬ 
ing will be required before any meaningful conclusions 
can be drawn. 

3.5. Modeling the mean spectrum 

The mean spectra shown in Figures [T] and [2] can be 
compared with earlie r UV spectra of NGC 5548 ob¬ 
taine d with HST fe.g.. lKorista et ai]ll995tlKaastra et al.l 
120141) . The earliest high-quality UV spectra of NGC 
5548 showed only weak absorption in the resonance lines 
(Lya, Nv, Siiv, and Civ), a factor that contributed to 
our selection of NGC 554 8 as a target for this investiga¬ 
tion. The 2013 spectra (jKaastra et al.l 120141) . however, 
revealed not only several strong narrow absorption 
features in the resonance lines, but also evidence for 
a relatively large “obscurer” that strongly absorbs the 
emission in the blue wings of the resonance lines. This 
feature is still present in our spectra obtained a year 
later, but is weaker than it was in 2013. The presence of 
this absorption, combined with the blending of various 
emission features (Lya and N v, Siiv and O iv], C iv and 
Hell), complicates analysis of these spectra. To charac¬ 
terize the emission-line structure of NCC 5548 and guide 
our selection of continuum windows for our emission-line 
flux measurements, we fit a heuristic model to the lines 


and continuum in the mean spectra. Our model for 
NGC 5548 is similar to that adopted by iKaastra et al.l 
(|2014D . and it includes the broad absorption features 
associated with all permitted transitions in the spec¬ 
trum. Our adopted continuum is a reddened power 
law of the form F\{X) = Fa(1000 A)(A/1000A)““. 
We correct for E(B — V) = O.OlTmag of Galac¬ 
tic extinction (iSc hl egel. Finkbeiner. fc DavisI 119981 
Schlufl^,^JjnLbehieijj2011[) using the prescription of 
Cardelli. Glavton. fc Mathia (119891) and Ry = 3.1. 
We do not apply any correction for possible inter¬ 
nal extinction in NGG5548. Longward of 1550 A, 
we also include blended Fe ii emiss ion as mod¬ 
eled by I Wills. Netzer. fc Wills! (I1985D . broadened 
with a Gaussian with full-width at half-maximum 
FWHM = 4000 km s~^. We model the emission lines 
with multiple Gaussian components. These are not an 
orthogonal set, and the decomposition is not rigorously 
unique, but they characterize each line profile well. 

For the brightest lines, we start with a narrow com¬ 
ponent, typically with FWHM rs 300 km s~^. This 
component is es sentially identica l to t he narrow com¬ 
ponent used by iCrenshaw et al.l (|2009t) for fitting the 
2004 STIS spectrum of NGG 5548. Since this narrow 
component is difficult to deblend from the broader com - 
ponents of each line, and since iGrenshaw et al.l (120091) 
saw little variation in narrow-line intensity over time, 
we fix the flux and widths of the narrow components of 
Lya, Nv, Giv, and Hell to the values we used to fit 
the 2004 STIS spectrum. The Siiv lines do not have a 
detectable narrow-lin e component. We note that while 
iPeterson et al.l (|2013D detected changes in the strength 
of the narrow [O ill] AA4959, 5007 lines over timescales of 
years, the variations over the last decade have been only 
at the few percent level. Similar variations in the narrow 
components of the UV lines would not be easily detected 
here because the narrow components are all so weak. 

Next we add an intermediate-width component with 
FWHM ss 800 km s~^ without ascribing physical mean¬ 
ing to it, using the STIS 2004 spectrum as a model. 
An intermediate-width component is included in Lya, 
Nv, Siiv, Civ, and Heii, as well as in the fainter 
lines C III* A1176, Si ii A1260, Si ii-kO i A1304, C ii A1335, 
N iv] A1486, O III] A1663, and N iii] A1750. For the weaker 
lines, this is often the only component detected, so its 
flux, width, and position are all allowed to vary. For the 
stronger lines, as for the NLR components, we again keep 
the fluxes and widths of the intermediate-width compo¬ 
nents fixed at the values found fo r the STIS 2004 spec¬ 
trum since ICrenshaw et al.l (|2009t ) found that these com¬ 
ponents vary only slightly in flux over several years. For 
the stronger lines, we next include broader components 
with FWHM ss 3000, 8000, and 15,000 km s~^, respec¬ 
tively. For the doublets of N v. Si iv, and C iv, we assume 
the line-emitting gas is optically thick and fix the flux 
ratio of each pair to 1:1, although for the 15,000 km s“^ 
component, only a single Gaussian is used. Finally, as 
can be seen in the RMS spectrum in Figure [21 there are 
two weak bumps that appear on the red and blue wings of 
the Giv emission-line profile at 1554A and ^ 1604A 
in the observed frame. These bumps are also present in 
the mean spectrum; we include a single Gaussian com¬ 
ponent to account for each of these bumps. 
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Figure 4. Interpolated cross-correlation functions (ICCFs) computed by cross-correlating each emission-line light curve with the 1367 A 
continuum light curve. The solid line is the ICCF for the entire data set. The dashed line is for the T1 subsample (the first 85 visits) and 
the dotted line is for the T2 subsample (the last 85 visits). 
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Figure 5. Velocity-binned results for the CivA1549 emission line. The top and middle panels show the mean and intrinsic RMS ((Jq) 
spectra, respectively. The bottom panel shows the centroid of the cross-correlation function for each individual velocity bin, corrected to 
the rest frame of NGC 5548. Velocity bins with no lag measurement contain too little total flux to reliably characterize the variations. In 
each case, the black line represents the entire dataset, and the T1 (first 85 visits) and T2 (last 85 visits) subsets are shown in gray and 
orange, respectively. Note that the shortest lags are found for the highest-velocity gas. 
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Figure 6. Velocity-binned results for the LyQ:A1215 emission line. The top and middle panels show the mean and intrinsic RMS (o-q) 
spectra, respectively. The bottom panel shows the centroid of the cross-correlation function for each individual velocity bin, corrected to 
the rest frame of NGC 5548. Velocity bins with no lag measurement contain too little total flux to reliably characterize the variations. In 
each case, the black line represents the entire dataset, and the T1 (first 85 visits) and T2 (last 85 visits) subsets are shown in gray and 
orange, respectively. The large gap at around —5000 km s“^ avoids the region of the spectrum affected by geocoronal emission and Galactic 
absorption. 


As described bv iKaastra et al.l (I2014li . we use an asym¬ 
metric Gaussian with negative flux to model the broad 
absorption troughs. The asymmetry in these Gaussian 
profiles is introduced by specifying a larger dispersion 
on the blue side of line center than on the red side. The 
asymmetry is fitted as a free parameter, and the resultant 
absorption line has a roughly rounded triangular shape 
with a blue wing extending from the deepe st point in the 
absor ption profile (for an illustration, see iKaastra et al.l 
l2014f) . During the first part of our reverberation cam¬ 
paign, when these absorption features were strongest, an 
additional depression appeared on the high-velocity blue 
tail of the main absorption trough. We use a single, 
symmetric Gaussian to model the shape of this addi¬ 


tional shallow depression. For absorption by N v, Si iv, 
and G iv, since the individual doublet profiles are unre¬ 
solved, we assume the lines are optically thick, so each 
line in the doublet has the same strength and profile. 

As a final component, we include absorption by 
damped Galactic Lya with a column density of V(H i) = 
1.45 X 10^°cm“^ (|Wakker. Lockman. fc BrownI l20ll]l . 
The full spectral model for NGC 5548, excluding the 
narrow absorption, is shown in the upper panels of Fig¬ 
ure [1] and Figure [21 superposed on the observed mean 
spectra 0 In future papers, we will apply this model to 

The complete model of G iv showing each of 
the individual components appe ars in the S upple men- 
tary Materials that accompany IKaastra et al.l II2Q14I) at 
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individual spectra to isolate the individual emission-line 
fluxes and to study absorption-line variability. 

4. DISCUSSION 

To put the results reported here in context, we note 
that NGC5548 has been monitored in the UV for RM 
pur poses on two previous occasions, as noted in Section 
Ol In 1989, NGC5548 was ob served once every four 
days for eight months with lUE (|Clavel et al.lll99lT ). In 
early 1993, it was observed every other day with lUE 
for a period of two months, and during the latter part 
of that campaign, it was also observed dai ly for 39 days 
with the HST Faint Object Spectrograph (iKorista et al.l 
Il995f) . The primary goals of these two experiments were 
quite different: the 1989 campaign was the first mas¬ 
sive coordinated RM experiment and it was designed to 
measure the mean lags for the strong UV lines. The 
1993 campaign was a higher time-resolution experiment 
that was designed to eliminate ambiguities from the 1989 
campaign. Specifically, its goals were: 

1. To measure the lag of the most rapidly responding 
line. He ii A1640. 

2. To determine whether or not there is a lag between 
the UV and optical continuum variations. 

3. To determine whether the wings and core of the 
G IV emission line have different lags. 

The first of these goals was met, but in the case of the 
other two, the data only hinted at results that are being 
confirmed by this project (Paper II and Figure [SJ. 

In addition to these RM programs, several HST GOS 
spectra of NGC 5548 were obtained in 2013 with the pri¬ 
mary goal of studying absorption features in the UV as 
support for an intensive X-ra y monitoring progra m un¬ 
dertaken with XMM-Newton (iKaastra et alJ[2r)i4[) . Our 
own results on variable absorption features constitute an 
extension of that effort and will be the subject of a future 
paper. 

Again, for broader context, during the AGN STORM 
campaign NGC 5548 was at about the same mean con¬ 
tinuum luminosity as it was during the 1989 campaign 
(but with a somewhat lower amplitude of variability), 
somewhat brighter than in the 1993 campaign, and de¬ 
cidedly brighter than it was in 2013, which was near the 
end of a lower-than-n ormal state that lasted several years 
(jPeterson et al.l[2013fl . The resonance lines showed much 
more self-absorption in this campaign than in either the 
1989 or 1993 observations, but less than seen in 2013. 
The emission-line lags were somewhat larger during the 
1989 campaign and the emission-line fluxes were higher, 
at least in part on account of much lower absorption in 
1989. The 1993 emission-line lags were similar to those 
obtained in this experiment, but again the line fluxes 
were larger, but less self-absorbed. 

As already noted, the response of the emission lines be¬ 
comes complicated during the second half of the present 
campaign. The He ii A1640 light curve seems to match 
the 1367 A continuum most closely; this line responds 
primarily to continuum emission at A < 228 A, imply¬ 
ing that the variations in the 1367 A continuum provide 

http: / / www.sciencemag.org/cgi / content / full / science. 1253787/DCl 


a reasonable proxy for the behavior of the hydrogen- 
ionizing continuum (A < 912A). Hell arises closer to 
the central source than the other emission lines, and it 
is also the only non-resonance line. More detailed analy¬ 
sis will be undertaken once the HeiiA4686 and Balmer- 
line results become available from our contempora neous 
ground-based monitoring program (jPei et al.ll2015li . 

In addition to determining the geometry and kinemat¬ 
ics of the BLR, we also wish to use these data to improve 
on previous estimates of the mass of the central black 
hole. However, the strong absorption in the blue wings of 
the resonance lines, which was very weak if even present 
in the 1989 and 1993 campaigns, preclud es using simple 
meas urements of the RMS spectra (e.s.. iPeterson et ^ 
l2004|i to make a mass estimate. More detailed modeling 
that we hope will lead to a more accurate black hole mass 
will be undertaken in future papers. 

To summarize briefly, we have presented the first re¬ 
sults from a UV spectroscopic RM study of NGC 5548 
undertaken with HST COS in 2014. We detect strong 
variations in the continuum and find clear delayed re¬ 
sponse of the strong emission lines, Lya, N v, Siiv, C iv, 
and He ii. A preliminary investigation shows that there 
is indeed a strong velocity-dependence of the emission¬ 
line lags, at least in the case of Lya and C iv, although 
blending and strong resonance absorption will make in¬ 
terpretation challenging. However, we have also shown 
that a heuristic multicomponent model can account for 
virtually all the spectral features. In future contribu¬ 
tions, we will use this model as a starting point to ex¬ 
plore the UV spectral variations in detail. We will also 
undertake a similar analysis of contemporaneous optical 
spectra in an effort to more completely understand the 
BLR geometry and kinematics. In the accompanying Pa¬ 
per II, we combine the continuum light curve presented 
here with high-cadence observations with Swift for a sim¬ 
ilar reverberation study of the accretion disk structure in 
NGG5548. 
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